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A relation between equilibrium, steady-state, and waiting-time dependent dynamical two-time 
correlation functions in dense glass-forming liquids subject to homogeneous steady shear flow is 
discussed. The systems under study show pronounced shear thinning, i.e., a significant speedup 
in their steady-state slow relaxation as compared to equilibrium. An approximate relation that 
recovers the exact limit for small waiting times is derived following the integration through tran- 
sients (ITT) approach for the nonequilibrium Smoluchowski dynamics, and is exemplified within 
a schematic model in the framework of the mode-coupling theory of the glass transition (MCT). 
Computer simulation results for the tagged-particle density correlation functions corresponding to 
wave vectors in the shear-gradient directions from both event-driven stochastic dynamics of a two- 
dimensional hard-disk system and from previously published Newtonian-dynamics simulations of a 
three-dimensional soft-sphere mixture are analyzed and compared with the predictions of the ITT- 
based approximation. Good qualitative and semi-quantitative agreement is found. Furthermore, 
for short waiting times, the theoretical description of the waiting time dependence shows excellent 
quantitative agreement to the simulations. This confirms the accuracy of the central approximation 
used earlier to derive fluctuation dissipation ratios (Phys. Rev. Lett. 102, 135701). For interme- 
diate waiting times, the correlation functions decay faster at long times than the stationary ones. 
This behavior is predicted by our theory and observed in simulations. 

PACS numbers: 82.70.Dd, 64.70.P-, 05.70.Ln, 83.60.Df 



I. INTRODUCTION 

The application of shear flow to dense liquids can 
dramatically change their transport and relaxation pro- 
cesses. Even if the timescale set by the applied shear rate 
7 is slow compared to a typical single-particle relaxation 
time To (i.e., the Peclet number Peo = 7T0 <C 1), it can 
interfere with and supersede the slow relaxation times T a 
of the system (i.e., the "dressed" Peclet, or Weissenberg 
number Pe = ^T a 3> 1). In this case, the slow relax- 
ation of the system is usually found to be accelerated 
by the shear flow, a phenomenon known from colloidal 
suspensions as shear thinning, because a pronounced de- 
crease in the apparent viscosity results. To this decrease 
corresponds an increase in the single-particle diffusivities. 
These changes in transport processes are found even if the 
average static structure of the system (at least as mea- 
sured through two-point correlation functions) changes 
only slightly. 

The drastic change from equilibrium to steady-state 
transport properties begs the question about transient 
dynamics: what happens if such a slowly relaxing liquid 
is suddenly subjected to shear, regarding its dynamical 
correlations as the system progresses from equilibrium to- 
wards its steady state? One can investigate these effects 
most easily by looking at the waiting-time dependent dy- 
namical two-time correlation functions: switching on the 
external shear flow at t = 0, one measures the correla- 
tions of dynamical variables between some waiting time 
t w > 0, and a correlation time t = t w + r > t w . Of 
particular interest are the so-called transient correlation 



functions, obtained for t w = 0, to be compared with the 
reference cases of equilibrium and steady state. 

In the present contribution, we address this issue by 
presenting an approximate relationship between these 
three relevant types of dynamical two-point correlation 
functions. 

Recently, the dynamical evolution after switching on 
shear flow of steady rate 7 has been addressed [H by a 
combination of techniques: theoretically, in the frame- 
work of mode-coupling theory (MCT) for colloidal rhe- 
ology and an integration-through transients (ITT) ap- 
proach, experimentally, using confocal microscopy, and 
with computer simulation for a damped Newtonian- 
dynamics model. All three methods yield a consistent 
picture: a shear stress a(t) builds up at t > after 
switch-on, but docs not grow monotonically towards its 
steady-state value CToo . It rather exhibits an intermediate 
"overshoot" at times corresponding to an overall strain 
7 = 7< w 0.1. Such stress overshoot phenomena are in 
fact well-known 0-111, despite their ubiquity, their 
microscopic origin, in particular for the fully homoge- 
neous flow profiles studied, remains somewhat vague. In 
Ref. simulations were able to connect it to a sud- 
den change in the mean-squared displacement (MSD) 
of a tracer particle: even for the directions perpendic- 
ular to the shear flow where no explicit advection oc- 
curs, one observes a super-diffusive regime as the tran- 
sient MSD, Sr 2 (t,t w = 0), leaves the equilibrium curve 
around 7 rj 0.1, to cross over to the (much larger) 
steady-state curve which it reaches at 7 « 1. In this 
regime, MD simulations found motion to be almost bal- 
listic, 8r 2 (t w OT/7,0) w t x with x ~ 2. In experiment 
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(closer to Brownian dynamics), this supcrdiffusion was 
not as pronounced, yielding x only slightly larger than 
1. Within MCT-ITT and an additional ad-hoc approxi- 
mation akin to a generalized Stokes-Einstein relation, it 
could be shown that the stress overshoot and supcrdiffu- 
sion are directly connected and originate from an over- 
relation of microscopic stresses: the transient stress au- 
tocorrelation function (called a dynamical shear modu- 
lus) does not decay monotonically to zero, but exhibits 
a "dip" in the corresponding strain regime where it be- 
comes slightly negative just before reaching its zero long- 
time limit. 

Thus, the details on the evolution from equilibrium dy- 
namics to far-from-equilibrium dynamics under shear are 
encoded in transient correlation functions. This raises 
a two-fold interest in these transient correlation func- 
tions: first, given information on both equilibrium and 
the steady state, what can one infer about the transient 
dynamics? Second, recalling that MCT and ITT generi- 
cally build upon the transient correlation functions, can 
one test their generic implications? 

The relationship among the various two-point dynam- 
ical correlation functions that we present in the follow- 
ing, builds upon the ITT formalism, without explicit ref- 
erence to MCT. It should thus hold quite generally, at 
least qualitatively. We demonstrate this by comparing 
with computer-simulation data for both a Brownian and 
a non-Brownian system. For dense liquids and colloidal 
suspensions in equilibrium, it is a well-tested paradigm, 
that the long-time behavior of the correlation functions 
does not (up to an overall time unit) depend on the type 
of short-time motion, be it ballistic (Newtonian dynam- 
ics) or diffusive (Brownian colloidal particles). This in 
fact defines the regime of quiescent "structural relax- 
ation" , where slow relaxation processes arising from col- 
lective caging of particles govern the dynamics of the sys- 
tem. 

This equivalence of Newtonian and Brownian systems 
does not need to hold far from equilibrium, although 
similar shear-thinning effects are seen both in colloidal 
suspensions and atomistic metallic melts. Indeed, differ- 
ences were observed regarding the extent of supcrdiffusivc 
motion in the transient MSD Note that in the MD 
simulations of Ref. [![, shear was implemented through 
the boundaries of the (periodically repeated) simulation 
box only, by Lees-Edwards boundary conditions. Hence, 
particles at the center of the box remain at rest also 
for a short time after switching on the flow, until a lin- 
ear shear profile propagates from the boundaries towards 
the center. It was argued that this time scale is short 
compared to the jt ~ 0.1 of interest. In Brownian dy- 
namics, the issue can be set aside, as there one immedi- 
ately modifies the solvent flow profile throughout the box 
when implementing shear (in addition to Lees-Edwards 
boundary conditions). This situation also being closer to 
what MCT-ITT models, is another motivation to comple- 
ment the MD simulation data already partially discussed 
in Ref. [l[ with new simulations incorporating diffusive 



short-time motion. 

The peculiar features of the transient correlation func- 
tions not only highlight possible differences in the non- 
equilibrium response of the different system types. They 
also provide direct tests of the MCT-ITT formalism. At 
the core of ITT is a reformulation of non-equilibrium av- 
erages in terms of history integrals over equilibrium av- 
erages where the full nonequilibrium time-dependence is 
kept in the evolution of the dynamical variables. In the 
case of steady shear, this is precisely the transient corre- 
lation function measurable in experiment or simulation 
where a shear flow is switched on immediately at t = 0. 

While the transient (t w = 0) correlation function is the 
natural object to be treated in MCT-ITT, steady-state 
and more generally i^-dependent correlation functions 
have, in this approach, to be calculated afterwards. Our 
main result here enables MCT-ITT to do just that: by 
obtaining equilibrium and transient correlation functions 
from any theory, correlators for all t w can be expressed. 

The paper is structured as follows: after some nota- 
tional clarification (Sec. [TTJ) , wc present in Sec. IIHI the 
main theoretical derivation of our formulas. Section IIVI 
is devoted to a schematic-MCT illustration of the results, 
while Scc.|V]presents the analysis of computer-simulation 
data. Section [VTI concludes the discussion. 

II. CORRELATION FUNCTIONS 

Given two dynamical fluctuating variables Sf and 5g, 
i.e., functions with zero average that depend on the state- 
point r of the system, one defines the two-point correla- 
tion function C fa (t, t w ) for t > t w 

C^(t,t w ) = Jj dTdT'Sg{T)P{T,t\T',t w )x 

5r(T')* tw (T'). (1) 

Here, P(T,t\T' ,t w ) denotes the conditional probability 
that the system resides at state point T at time t = t w +r, 
given it was at state point V at time t w . ^> tro (r') is the 
probability that the system is at V at time t w , and of 
course just the (nonequilibrium) distribution function. 
The latter is assumed to be equal to the equilibrium 
distribution for times t < 0, v E't<o(r) = ^e(r), and to 
asymptotically reach a time-independent steady state for 
long times, ^-^(IT) = ^ s (r). These reference limits de- 
fine the equilibrium and steady-state ensemble averages, 

(■■•>= / dT* e (T)... , (2a) 

(•••) w = / dr*,(T)... . (2b) 

The conditional probability P in Eq. ([I]) encodes the 
dynamics of the system, subject to the external field. We 
will use the following property Q : for the case of switch- 
ing on a constant shear flow, P(T, t w + t\T', t w ) becomes 
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FIG. 1: Schematic representation of the times appearing in 
the two-point correlation functions, Eq. {T}: waiting time 
t w > 0, measurement time t, and correlation time r = t—t w > 
0. t = corresponds to the time where homogeneous, linear 
shear flow is instantaneously switched on. 



independent on t w for all t w > 0. For Brownian dy- 
namics, this means that both ^t w and P obey the same 
differential (Smoluchowski) equation, i.e., the underly- 
ing stochastic process is assumed to have a Markovian 
property [f|. As long as ^ tm (T') evolves, the correla- 
tion function C^ 9 (t,t w ) will depend on its two time ar- 
guments separately; we will generally call these functions 
"waiting-time dependent" (for waiting time t w ). When it 
is clear from the context, we omit the superscript denot- 
ing the variables and abbreviate Ct w (t) = C(t w + r,t w ). 
While potentially interesting, we ignore correlation func- 
tions formed with t w < and t > 0. Figure [1] schemat- 
ically summarizes the sequence of correlation and mea- 
surement times. 

As t w — > oo, the steady-state correlation function is ap- 
proached, which we denote by Coo(t) = C(t, +oo). Cor- 
respondingly, as t w < —t, one obtains the equilibrium 
correlation function, C e (r) = C(t, — oo). Both Coo(t) 
and C e (r) are functions of the time difference r only. 
Among the general waiting-time dependent correlation 
functions, a particular role is played by the t w = case: 
recalling = we recognize the so-called transient 
correlation function Cq(t) = C(t, 0), where the time evo- 
lution as determined by the transition rates P is the 
noncquilibrium one, but averaging is performed with the 
equilibrium distribution. 

Note that we do not in general assume C e (r) to decay 
to zero for r ^ oo: in (idealized) glass states, it attains 
a finite positive long-time limit, called the noncrgodicity 
factor or glass form factor Q. 

In comparing with computer simulations, we set for 
simplicity / = g to be the one-particle microscopic num- 
ber density, 5f = cxp[iq • r s ] with wave vector q and 
the position of the singlcd-out particle r,,. To further 
simplify the discussion, we restrict ourselves to wave vec- 
tors perpendicular to the flow direction. This obliterates 
the need of introducing wave- vector advection in order to 
account for the affine motion imposed by the shear Q . 

Connected to the zero-wavevector limit of tagged- 
particle density fluctuations is of course the mean- 



squared displacement (MSD), 

6r 2 (t,t w ) = Jl dTdT' [r s (r)-r s (r')] 2 x 

P(T,t\T',t w )* tw (T'), (3) 

also schematically written as Sr 2 (t,t w ) = ([r s (t) — 
r s(t w )] 2 ) sl where (• • denotes averaging over particles 
and runs as done in a simulation. Again, the MSD comes 
in its equilibrium (5r 2 (T)), steady-state (^^(t)), and 
transient (Sr 2 (r)) varieties. 

We generally consider a system of TV spherical particles 
without internal degrees of freedom, enclosed in a volume 
V. Choosing units, we set the thermal energy k^T = 1 
throughout. For the stochastic-dynamics simulations as 
well as the theory, we assume diffusive short-time motion 
governed by a bare diffusion coefficient Dq = 1; we do 
not take into account explicit solvent or hydrodynamic 
interactions among the particles. A typical interaction 
diameter of the particles, a = 1, sets the unit length. 

We choose coordinates such that the external flow acts 
in the Cartesian x-direction (called flow direction), and 
varies along y (called gradient direction). In 3D, the 
system is invariant along the z-axis (neutral direction). 
Thus, the velocity field induced by the shear can be 
written as v(r) = •jySc with the velocity-gradient ten- 
sor k = 7xy, where x is a unit vector in the direction of 
x. 



III. MICROSCOPIC THEORY 

Within the ITT formalism, we now derive expressions 
for the time-dependent correlation functions of interest. 
We begin by recalling the exact starting points of ITT, 
before introducing approximations that lead to our final 
result, presented in Sec. IIII CI 

A. Integration Through Transients 

For the theoretical derivation, assume interaction 
forces among the particles to be F t = —diU(i = 
1,...N), where U is the total potential energy of the 
system. In the thermodynamic limit, the particle distri- 
bution function ^(r) of the Brownian system subject 
to homogeneous shear flow described by the velocity- 
gradient tensor K:(t)j_is then taken to obey the Smolu- 
chowski equation [9|, [HJ , 

ae* t (r) = n(t)* t (r), (4a) 
n(t) = n e + ssi(t) = J2 9 >- & - F > - K ^ ■ r *l • ( 4b ) 

i 

Here, ft is the Smoluchowski operator (SO), consist- 
ing of the equilibrium (quiescent) contribution, f2 e = 
di ■ [di — Fj] , and the noncquilibrium term represent- 
ing homogeneous driving. For the case considered here 
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(switching on constant shear flow of rate 7 at t = 0), 
<5fi(i) = 50 = — J+ di ■ k ■ Yi independent on t for t > 
and zero else. Hence, 



O e for t < 0, 
fi w for i > 0, 



(5) 



where does not depend on time. The equilib- 

rium distribution function "J e is the stationary solu- 
tion of Eq. (j4]) without shear, D, e *f> e = 0, viz. ^ e oc 
cxp(— U/ksT). Including shear, ty s = defines the 
steady-state distribution. In this stationary state, the 
distribution function is time-shift invariant as in equilib- 
rium, but the system is not in thermal equilibrium due 
to a non- vanishing probability current (Tlj . 

The integration through transients (ITT) formalism al- 
lows to reformulate the nonequilibrium averages formed 
with the (unknown) t^-depcndcnt distribution $ tii in 
terms of equilibrium averages. Formally solving the 
Smoluchowski equation, Eq. (j4|), as an integral equation 
with the boundary condition ^(t = 0) = fy e , one gets for 
t w > 0, 



ds fi™ 



„fi«> : 



*e- (6) 



Recalling * e = SQty e = a xy ^ e [l(J, integration by 
parts yields 



dT^ t 



dTV e 



1+7 



ds a xy e 



n 1 ' 1 



(7) 



with the microscopic (potential) stress tensor element 
<r*y = -J2iFiVi M- Here, ft* = J2 i [9 i +'F i +v i -K T }-d i 
is the operator adjoint to fi™. With t w — > 00, the 
steady-state average in Eq. ^ immediately follows 



v(7) 



7 / ds ((7rr,„e 



(8) 



The Smoluchowski equation, Eq. Qj, is also taken 
to determine the conditional probability P appearing in 
Eq. CQ), and thus (EGJJ, for t w > 0, P(T, t w +t\V, t w ) = 
exp[0 W r]5(r - T'). Rewriting in terms of the adjoint 
operator, this gives 



Cf 9 (T) = (6f*e Q ' T 6g) 



(7) 



(9a) 



and, for the transient correlation function, 

C /ff (r) = (tffe^) . (9b) 

Note its distinction from the equilibrium correlation func- 
tion, where the equilibrium adjoint SO, Q\ 7 appears, 



(9c) 



For the general two-time correlation function at finite 
t w > 0, inserting into Eq. ([7]) gives 



7 J ds (a xy e' l ' s 5te nW 5 9 y (10) 

This equation is formally exact, although the evaluation 
of the dynamical three-point average in the integral will 
generally be hard. We are therefore forced to introduce 
approximations at this point. 



B. Approximations for Correlation Functions 

To simplify the discussion, we now restrict ourselves to 
auto-correlation functions (Sf = Sg) of dynamical fluctu- 
ations without explicit shear-advection, Sf = f({yi, z.i}). 
Similar results can be expected for correlation functions 
involving shear-advected quantities (i.e., dynamical vari- 
ables explicitly depending also on positions Xi along the 
shear flow), but one then has to be careful in first ex- 
tracting the affinc transformations induced by the steady 
shear. 

To obtain a tractable expression for the general in- 
dependent correlation function, let us apply a famil- 
iar identity in the Zwanzig-Mori operator formalism (cf. 
Eq. (11) in Ref. [TJ| and also Ref. [UJ]): introducing a 
projector onto Sf, V f = Sf)(Sf*Sf)~ 1 (Sf*, with com- 
plement Qf = 1 — Vf, we get from Eq. (fTOjl. 



C/»=Qf(T) 



1 + 7 



tw 



(Is 



d O t s 



Sf*Sf 



(Sf*Sf) 



7 / dr 1 

>o Jo 



ds 



(7 XV C 



W s 



Sf*U(r 



(Sf*Sf) 



')Sf) 



(11) 

with the restricted time evolution operator 

U{a) = Q f exp[Q f rtQ f a}Q f rt . (12) 

We thus identify two contributions to the difference 
between the nonequilibrium waiting-time dependent cor- 
relator and the transient one. The first is a r-independent 
renormalization of the equal-time value and corresponds, 
e.g., to the difference of distorted and equilibrium static 
structure factor if t w — > 00 0, [JJ[ ■ Note that in Ref. [1] , 
only this term for the difference of the correlators is con- 
sidered. It vanishes for tagged-particle density fluctua- 
tions, since with Sf = exp[iq • r s ], the average of Sf*Sf 
is unity in any ensemble. 

The second term contains a more complicated depen- 
dence on both t w and r, and cannot easily be evaluated. 
But one recognizes that inserting a projector onto a xy 
before the Sf* term allows to factorize the integral ac- 
cording to the different time dependences. The right- 
hand part containing U(t — r') then becomes propor- 
tional to (<r xy Sf*U(T — r')Sf), for which the operator 
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identity that led to Eq. (fTTj) can be rolled back by not- 
ing {a xy Sf*Sf} = due to symmetry (5f*Sf is symmet- 
ric in coordinates x and y, while a xy is antisymmetric). 
Thus, assuming the dominant part of the last integral in 
Eq. (fTTj) to be given by the projection of Sf onto a xy , we 
get 



C{Jt) ps a s {t w )C J (T) + ja(t w ) (a xy Sf*e u ' T Sf 



<I 



(13) 



where we have abbreviated the static renormalization by 
a/(t„L= 1 + 7 Io W ds (a xv exp[n 1 s]6f*6f)/(Sf*Sf). In 
Refs. [HI, HBl the simplified version with af(t w ) = 1 was 
considered which, again, holds exactly for tagged particle 
dynamics. 

Note that formally, we had to introduce a projector 
V a = cT X y}(a xy a X y}~ 1 (a X y that is ill-defined for the case 
of hard spheres, as there, the instantaneous shear modu- 
lus (<7 xy G xy ) diverges 0, 0|. However, we only require 



a(t w ) 



{o~ X yO~ xy ) 



■ ds . 



(14) 



i.e., the integrated normalized shear modulus [10t ll8H21l ]. 
to exist. We assume that this integral can be regularized 
for hard spheres, as outlined in Appendix 1X1 

The remaining correlation function in Eq. (|13[) is noth- 
ing but the waiting-time derivative of C{ (r) at t w = 
[ll| . as is immediately clear from taking the ^-derivative 
on both sides of Eq. (TTT 



, xy 5fe n ^Sf) = 



_d_ 

dt„ 



(15) 



It describes the initial change of the two-time correlator 
with t w at fixed correlation-time window r. Our approx- 
imation then reads 



C' (t) ps a f (t w )Ci(r) + a(t w ) J-C£ (r) 



dt„ 



, (16) 



and Eq. (|13[) can be interpreted as "coupling at t w - 
0" , incorporating the exact result to first order in t n 
recalling a(t w ) = t w + 0{t 2 w ) and a(t w ) = 1 + G{t 2 w ), 



C/ iu (r)=C / (r)+A C /jr) 



t w +G{{it w f) . (17) 



Equation (|13l) extends this identity to finite t w by ac- 
counting for the static change otf(t w ) exactly (in princi- 
ple), and relating the further t w dependence to the inte- 
grated shear modulus a. 

We still have to close this approximation by relating 
the waiting-time derivative to known correlation func- 
tions. As was shown in Refs. [TJ IH|i one can j using 
integration by parts and the identity <5fi^ Sf = 0, arrive 



at 



Ot, 



Sf* S& e n+ T Sf 



dr 



Gt{r)-(srnle a ^Sf) . (18) 



This equation highlights the connection of the waiting- 
time derivative to time derivatives of correlation func- 
tions: the derivative of the transient correlator Cq(t) 
has two parts, one containing the equilibrium operator 
fig, and one containing the nonequilibrium shear- induced 
<5fi^. The former term corresponds to the short-time dy- 
namics of the correlation function, unaffected by shear as 
long as Peo -C 1, while the latter term, the waiting-time 
derivative, is governed by the shear-induced decay of the 
correlator at long times. 

The equilibrium derivative fig Sf* in the last term of 
Eq. (jTHJ) de-correlates quickly as the particles loose mem- 
ory of their initial motion even without shear. In this 
case, the latter term is the time derivative of the equilib- 
rium correlator, C{ (r). A shear flow switched on at r = 
will generally lead to even faster decorrelation, prompt- 
ing us to approximate e° + T ps e f2 = T Vf e~ n = T e at r . This 
approximation used in the last term in Eq. (|18[) as well 
as in Cq (r) leads to 

(Sf* fig explfit r]5f) (sf* fi e exp[fi e r]Sf) 

ps (19) 

(6f*exp[tfT]8f) (<J/*exp[fi e r]5/) 



and thus 



Sf* fij e ot T 5f 



Inserting this approximation in Eq. (|18|) yields [11| 



C[{t) dr 



_d_ 

dt„ 



C((r) 



d_ 



Cl(r) dr AT) 



(20) 



(21) 



The two terms in this equation have an intuitive inter- 
pretation: if 7T < 1, there holds Cjf (r) = C{ (t) + 0(7t) 
[17| , and the right hand side of Eq. (j2Tj) cancels in leading 
order in jr. This is expected on physical grounds, since 
the short-time decay of the correlation function is inde- 
pendent of t w at least for small t w . For r = 0, Eq. (|2Tj) 
yields zero exactly, in agreement with Eq. (fTS"]) . where 
(cr xy Sf*Sf) = due to symmetry. On the other hand, 
for 7T = 0(1) with Pe ^> 1 (i.e., the relaxation time of 
C e (r) is much larger than the shear- induced relaxation 
time 0(1/7)), the last term in Eq. (j2"T)) vanishes, and the 
waiting-time derivative is given by the time derivative of 
the transient correlator. We thus refer to the term on the 
left hand side and the last term in Eq. ([21]) as long-time 
and short-time derivatives, respectively. 

The approximation leading to Eq. ()21[) can further be 
made plausible by considering states that are glassy in 
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the quiescent equilibrium; setting the second term on the 
right-hand side to zero, and writing out the derivatives, 
one gets 



t — v oo . 



(22) 



for small St. This embodies the physical argument that at 
large times, whenever the equilibrium dynamics is frozen 
and the transient correlator is on the plateau, shear ef- 
fects set in as function of t rather than r. 



Relation for the Two-Time Correlator 



Equations (|21[) and (|16|) taken together yield an ap- 
proximation that allows us to study the waiting-time 
dependence of the non-equilibrium two-time correlation 
function, 



a f (t w ) 
d 



a(t w )-[ln\c£(T)\-ln\Cf(T)\ 



(23) 



A brief discussion of this result might be in order. First, 
we recognize that for weak shear, Pe < 1, the second 
term in Eq. (|23[) does not contribute, as in this regime 
Cq (t) » C/(r), and hence the normalized waiting-time 
dependent correlation function likewise does not change. 
This ensures that we correctly recover linear response. 
On the other hand, for Pe > 1, the derivative of the 
transient correlation function will dominate the second 
term in the equation. 

Let us also note the equivalent approximation to 
Eq. (|23|) for the mean-squared displacement, easily de- 
rived from the q — > limit of the corresponding tagged- 
particle density correlation function. Considering the y- 
or z-direction for simplicity, 



Sz (t) = lim 



l-C'(r) 



q-S-0 



(24) 



where 5f = cxp[iqz s ] with z s the ^-coordinate of the 
tagged particle. Performing the q — » limit in Eq. (f2"3")l 
directly yields (recall that aif(t w ) = 1 for tagged-particle 
density fluctuations) 



6zt(r) « 5zl{r) + a(t w )- (<5z 2 (r) - <5z e 2 (r)) . (25) 

Considering only directions perpendicular to the shear 
direction, the MSD is linear in time for long times, i.e., 
Sz 2 (t) ~ t as r — > oo. If 7T ^> 1 (and thus jt S> 1), it 
is plausible that the transient MSD describes the same 
diffusivity as the stationary one, and the time-derivatives 
of the two functions have to be equal in that limit, as 
reproduced by Eq. (f2"5| . 

Equations ([23")) and (f2"5|) constitute our main theoret- 
ical result: calculating, e.g., Cq(t) and C/(r) within 



MCT, the above equations give access to the general 
waiting-time dependent C{ (t), including the steady- 
state correlation function usually measured in experi- 
ments. 

The equations can also be rewritten in order to de- 
termine the transient correlation function from the more 
commonplace equilibrium and steady-state ones. To this 
end, note that Eq. (|2"3")l is solved by 



C£( T ) = e -«/(*« W*(*«)c£(t) 

Jo <n*u,) cl{s) 



(26) 



Taking t w — > oo yields an expression determining the 
transient correlator Cq(t) in terms of both the equi- 
librium one, C/(r), and the stationary one, C^,(r). 
For small waiting times, <r(t w ) = t w + 0(t^) can be 
used, and for tagged-particle correlation functions where 
otf{t w ) = 1, no unknown parameters remain in Eq. (|26p . 
For the mean-squared displacement, a similar transform 
holds. 



8zI(t)=8zI{t) 



,{s)-&zl{s) 



a{t w ) 



ds . (27) 



Incidentally, these forms are also, for direct testing with 
computer-simulation data, more stable numerically than 
Eqs. ([23jl and (|25|) as they do not involve time-derivatives. 
We will therefore use them primarily in Sec. IVl where we 
use C/(r) and C( (r) taken from computer-simulation 
data to assess the quality of the approximation by com- 
paring the calculated and simulated Cq (t). 



IV. A SCHEMATIC MCT MODEL 



Equations ([23]) and (f2"5|) represent the central theo- 
retical result of our paper. In this section, we visualize 
their physical content by choosing a simple toy model for 
calculating the transient two-time correlation functions. 
For the purpose of our discussion, it is easiest to choose a 
schematic model of mode-coupling theory, as these mod- 
els are very successful both in analyzing real-world data 
and for understanding the generic features of colloidal 
systems under shear. 



A. Schematic Equations 

Let us consider a single, normalized transient corre- 
lation function Cq(t) — <p(r), to represent the collec- 
tive density-fluctuation correlators for some dominant 
(nearest-neighbor) length scale. Recently, a schematic 
model allowing to treat arbitrarily time-dependent flow 
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has been proposed [22j , whose equations of motion in the 
case of steady simple shear reduce to 



= r o a T 0(r) + 0(t) + / dT'm(T,T-T')d T >4>(T')dT' , 
Jo 

(28a) 

m(a, 6) = hi(a)h 2 {b)m(b) , (28b) 
m(a) = [ Wl 0(a) +w 2 ^(a) 2 ] , (28c) 
fci, 2 (a) = 1/ [l + ( 7 a/7c) 2 ] . (28d) 

Here, the hi{a) are ad-hoc forms for the strain-induced 
reduction of the memory kernel, inspired by the appear- 
ance of similarly time dependent terms in the original 
MCT vertices due to wave-vector advection. Its precise 



form is not crucial, and we follow the choice of Ref. [22 1 
by choosing a simple decaying function that is even in 
the strain. There, j c = 0.1 was introduced to model 
the typical cage-breaking length scale: strains of about 
10% mark the point where noticeable strain reduction of 
memory effects sets in. 1/tq is an initial decay rate that 
serves to set the unit of time. 

For 7 = 0, Eqs. (f28|) reduce to the well-known F12 
model of quiescent MCT [23[, whose solutions provide 
the equilibrium correlator e (r). Setting h\ = 1 further 

reduces our model to the F^ model [l7^ . originally pro- 
posed for the analysis of steady shear flows. This model 
provides excellent fits to the flow curves from large scale 
simulations [24[, and the extension setting h\ = ft, 2 does 
not qualitatively change these flow curves, while keeping 
a closer connection to the more general time-dependent 

flows H3. 

The F12 model has glass transitions along a line of 
coupling parameters (v^v^), where the long time limit 
/ = lmv^oo 4> e (r) jumps discontinuously from zero to 
its critical value f c . The choice V2 = 2 is known to 
yield good agreement with the asymptotic features ex- 
pected for the hard-sphere glass transition; it implies 
vl = 2{\/2 — 1). The separation parameter e then 
serves to quantify the distance to the transition: we set 
(vi, V2) = (fi, + £ ) sucn tnat £ > indicates glassy 
states, e < fluid ones. 

In order to evaluate Eq. (f2"3"f. we further need a 
schematic- model version of Eq. (fFT| yielding a(t w ). We 

let unm 



<fi(s) ds . 



(29) 



which can be regarded as a schematic version of the gen- 
eralized Grccn-Kubo relation derived within ITT and the 
MCT approximation [22| , where we approximate the dy- 
namical shear modulus G(s) w 0(s)/3, a reasonable ap- 
proximation in particular at long times. The factor 1/3 
accounts for the fact that the plateau in the shear mod- 
ulus is empirically found to be smaller than that of the 
correlator. In fact, Eq. (f2T)|) neglects prefactors and an 
anisotropic wave-vector integral that appears in the mi- 
croscopic Green-Kubo relation; this means also that we 



have lost the correct description of the stress overshoot. 
To include this effect, one would need to evaluate Eq. (|2"9")l 
microscopically, as done e.g. in Ref. [l[. Eq. (|29|) is thus 
to be regarded merely as a plausible closure that incorpo- 
rates the structure of MCT-ITT that cr(t w ) is dominated 
by an integral over the density correlation functions, ef- 
fectively cut off by the slow relaxation time of those cor- 
relators. 

Equations analogous to Eqs. (|28|) hold for the 
schematic transient tagged-particle correlation function 
4> s (t); the only difference is in the precise form of the 
mode-coupling kernel m s (a), Eq. (|28c|) . Its microscopic 
expression for tagged-particle density fluctuations can be 
worked out [H], and will be discussed elsewhere. For our 
purpose, we copy the form of the well-known schematic 
quiescent tagged-particle model, the so-called Sjogren 
model, m s (a) = v s (f)(a)(j) s (a). Here, a coupling coeffi- 
cient v s > appears that describes the strength of the 
tagged-particle coupling to the collective density fluctu- 
ations. This parameter plays no qualitative role in the 
further discussion; we fix it to v s = 5. 

The transient MSD (in the neutral or gradient direc- 
tion) is the solution of a similar memory equation, cf. 

Refs. dm, 

t Sz 2 (t) + f m s (r - t')Sz*(t') dr' = 2r , (30) 
Jo 

where we denote the schematic-model transient MSD 
by Szq in order to avoid confusion with its microscopic 
counterpart. In principle, its memory kernel will not 
be identical to the one appearing in the taggcd-particlc- 
corrclator equation, but the MCT approximation for the 
self-density fluctuations imply that both these memory 
kernels are bilinear functionals of </> and s , so that they 
can be approximated as equal on the schematic level. 
Note however that rh s must be a single-time function, 
as is found in the full microscopic derivation 25|, in 
order to recover long-time diffusion. We therefore set 
m s (a) = h2{a)rh s (a) . 



B. Results of the Schematic Model 

We now turn to a discussion of the central Eqs. 
and (|25p with the aid of the schematic transient and equi- 
librium correlators and MSD defined in the previous sec- 
tion. In order to highlight the non-trivial effect of the 
waiting time, we set a = 1 now, keeping all correlators 
normalized to unity at r = t — t w = 0. Figure [5] shows 
results for (j> tm (r) , the schematic waiting-time dependent 
correlator calculated via Eq. (|23|) (replacing the general 
correlation functions C e and Ct w with the schematic ones, 
(j) e and 4>t w , omitting the ^-subscript only for the tran- 
sient function cf)). A glassy state was chosen, e = 10~ 3 , 
so that 4> e {T) attains a finite long-time limit /, and the 
decay of cj>t w (r) as r — >• 00 is solely due to the shear 
flow. In this case, the last term in Eq. (1231) vanishes 
at long times, and the time derivative of the remaining 
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FIG. 2: Dynamical two-time correlation functions, Eq. (|23|l . 
using the F^lj' model P3] as input for the transient correlator, 
in the glass (e = 10 -3 , see text) at shear rate jro = 10 -9 , 
for various waiting times t w as indicated (curves from top 
to bottom); thick lines indicate the transient (t w = 0) and 
steady-state (t w = oo) correlators. The dotted line represents 
the equilibrium correlation function. The lower panel displays 
the final decay as a function of strain -y(f — t w ). 



term is negative for all t. The decay of the transient 
correlator from the plateau can be well approximated by 
0(t) w /exp[— {^t/^cY] with exponent \i = 1.2, i.e. it 
shows "compressed exponential" behavior as a signature 
of the non-steady dynamics. Recall that in equilibrium 
colloidal suspensions, \x > 1 is excluded by the properties 
of the Smoluchowski operator [2?], HU. Approximating 
the decay for argument's sake as a simple exponential, 
we immediately see that 4't w ( T ) ~ <X T )(1 — j&itw)) for 
7T = 0(1). Since ct(oo) remains finite, Eq. (|23p indeed 
describes the asymptotic approach to a steady-state cor- 
relator (fioofr) as t w oo. As seen in the figure, this 
approach occurs on a time scale jt w 5%. The differ- 
ence between steady-state and transient correlation func- 
tions becomes noticeable only once the correlation func- 
tions decay from their plateau, at jr « 0.001; it vanishes 
as the functions decay to zero and is most pronounced 
around 7T w 0.01. 

In the liquid, e < 0, similar effects as those described 
above are seen, but only in a regime where Pe = 7T Q 1, 



FIG. 3: Transient (solid lines) and stationary (dashed lines) 
correlation functions of the schematic model for three state 
points: top curves demonstrate e = — 10 -3 (liquid with a 
relaxation time T a /ro = 0(1O 7 ), shifted by 0.8 vertically), 
middle curves e = (transition point, shifted by 0.4), bottom 
curves e = 10~ 3 (glass). Shear rates are Peo = 7T0 = lO^ 2 ™" 1 
with n = 1, ... 4 (left to right). 



where r a is a time scale characterizing the slow relaxation 
of the equilibrium correlator <fi e (r) (the "a" time scale in 
glassy liquids), defining the dressed Peclet number. In 
this case, again, the last term in Eq. (|23[) can be dropped 
as a small correction, as e (r) ps /, its plateau value, 
for r ?s I/7. This is visualized in Fig. [3l where results 
for various e and different shear rates 7 are shown. For 
Pe S> 1 (and Peo < 1), curves on both sides of the glass 
transition show qualitatively identical evolution with t w . 
However, in the liquid, a linear-response regime exists for 
Pe <C 1, where Eqs. (|2"8"|) describe a transient correlator 
that is itself unaffected by shear, <f>{r) ~ 4> e (i~). This is 
exemplified by the top right curve of Fig. [3J where a liquid 
state of the F12 model with small shear rate correspond- 
ing to Pe = ■yTa = 10~ 2 was chosen. Equation (f2"3")l then 
correctly describes the fact that all i^-dependent corre- 
lation functions are equal. Note that this holds possibly 
only up to a normalization expressed by a(t w ), reflect- 
ing the fact that the static structure may be distorted by 
shear. 

In the ideal glass, the quiescent relaxation time is in- 
finite, and hence no linear-response regime exists. In- 
stead, the transient correlator 4>( T ) decays on a time scale 
0(1/7) for arbitrarily small shear rates. Hence, for large 
enough jr and 7 — > 0, Eq. (|23[) becomes invariant under 
the transformation r 1— > jr and t w t— > jt w in the glass. 
This yields a nontrivial prediction, namely that the dif- 
ference between transient and steady-state correlators in 
the shear-molten glass does not vanish as 7 — > 0. Rather, 
a scaling limit is exhibited if one considers the correlation 
functions on rescaled times 77-, where both the transient 
and the stationary correlators attain (different) master 
curves. The approach to this scaling is demonstrated 
by Fig. [U where the e > curves of Fig. [3] are repro- 
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FIG. 4: Transient (solid) and stationary (dashed lines) corre- 
lators of the schematic model in the glass, e = 10 -3 , for shear 
rates Pen = lCT 2 ™" 1 , n = 1, . . . 4. 
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FIG. 5: Transient (solid) and stationary (dashed) mean- 
squared displacements in the schematic model, Eq. (|30[1 and 
Eq. (25) , for a glassy state, e = 1CP 3 , with shear rates 
Pen = 7T = 10 -2n_1 , n = 1, ...4 (from left to right). The 
quiescent equilibrium MSD is shown as a dotted line. 



duced as functions of jr. One clearly identifies the ap- 
proach to two master functions looking at the lowest two 
shear rates, jtq = 10~ 9 and jro = 10~ 7 . One notes 
also that the scaling is only approached for rather small 
shear rates, and thus not easily verified. Although this 
nonanalytic limit presents an interesting test of MCT, it 
will in general be difficult to probe in simulation, as one 
would need to construct transient correlation functions in 
the properly ensemble-averaged glassy initial state. The 
physical reason behind this nonanalytic limit is the exis- 
tence of a dynamical yield stress predicted by the theory, 
rendering the 7 — > limit under steady shear a singular 
one. 

Figure [4] demonstrates another prcasymptotic effect of 
the limit of vanishing shear rate: identifying the regime 
of accumulated strain where transient effects are largest, 
we found 7T w 0.01 from the small-shear-rate regime dis- 
cussed above. However, at larger shear rates, this regime 
shifts to almost jr w 0.1; these arc the 10% strain cor- 
responding to a typical localization length of hard-core 
particles in the glass, argued for in Ref. [lj. 



For the mean-squared displacements, via Eq. (|25|) . a 
very similar discussion holds, as exemplified by Fig. [5] 
The qualitative differences between equilibrium, tran- 
sient, and steady-state MSD curves is the same as for 
the correlators (although now, for obvious reasons, the 
functions increase, rather than decrease, with increasing 
t w for fixed r). The behavior is also in qualitative agree- 
ment with experimental and computer-simulation results 
of Ref. [l[. However, there, an intriguing super-diffusive 
regime for the transient MSD was observed which is miss- 
ing in our schematic model. It can be argued that this is 
due to Eq. (|30p , where we model the memory kernel of the 
MSD by a strictly positive function rh s (a), whereas the 
true memory kernel leading to a superdiffusivc regime 
should exhibit a small time window of negative values. 
Indeed, superdiffusivc behavior was found in Ref. [l[ to 
be connected with the stress overshoot phenomenon (not 
modeled in our schematic approach) that was argued to 
enter fh s {a) via a generalized Stokcs-Einstein approxima- 
tion. Thus, the qualitative agreement of the ^-evolution 
of MSD curves between our model and the simulation 
data highlights that the superdiffusivc motion is not nec- 
essarily connected to the physics of crossing over from 
equilibrium to steady state after switching on shear flow. 

V. COMPARISON WITH COMPUTER 
SIMULATION 

We now turn to a discussion of the waiting time de- 
pendence of Ct m (t) as found in computer simulation. 

We first report the findings from a stochastic-dynamics 
computer simulation for a two-dimensional system of 
hard disks. To avoid crystallization, an equimolar binary 
mixture with diameters a — 1 (taken as the unit length) 
and 1.4cr is chosen; this is the same system as studied 
earlier in steady state (29j |. The number density is con- 
veniently expressed as the packing (area) fraction of the 
system ip. Previous simulations found a glass transition 
at packing fractions around tp as 0.8. 

The simulation is modeled after the so-called event- 
driven Brownian-dynamics (ED-BD) algorithm [3(| ; this 
algorithm provides an approximate solution to the 
stochastic differential equation underlying Eq. (|4]) with- 
out shear flow, 

= -Qdx + dF r , (31) 

where x is the configuration-space vector describing the 
particle positions, ( is a friction coefficient and F T 
is a random white-noise force obeying the fluctuation- 
dissipation theorem, (F? (t)F* (f)) = 2k B T(5(t - t')Sij. 
Setting the amplitude of the noise correlation to unity 
fixes the unit of time. The interaction among the hard- 
sphere particles translates into boundary conditions that 
no two spheres overlap at any point in time. The ED-BD 
algorithm is, in essence, a rejection-free hybrid Monte- 
Carlo scheme that works by selecting a small time step 
At, during which a free Brownian particle undergoes 
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a displacement with variance (Ax 2 ) = 2DoAt in each 
Cartesian direction. Trial moves are first drawn accord- 
ing to the Gaussian distribution of that variance, and 
then corrected for unphysical overlaps. The overlap re- 
moval is performed by assigning to each originally drawn 
displacement Ax a tangent vector u = Ax/ At and a lin- 
ear curve parameter s £ [t,t + At]. Particles are then 
displaced along the tangent vectors by a linear mapping 
from s = t to s = t + At. Whenever at some s c two 
particles i and j start to overlap, the corresponding trial 
vectors Ui and Uj are reflected along the plane perpen- 
dicular to the particles' separation vector, which ensures 
no-flux boundary conditions on the spheres' surfaces if 
At is small enough. Effectively, this translates into per- 
forming "elastic collisions" with the Ui and Uj treated as 
velocity vectors. The procedure is continued, taking care 
of all of the possibly many s c in the same fashion until 
s = t + At. One is then guaranteed to have a new config- 
uration that is overlap-free and that the phase space is 
sampled ergodically by diffusive motion of a free diffusion 
coefficient D = At/2. [3rj| . 

Linear shear flow is incorporated in this algorithm fol- 
lowing Ref. [29[ by shifting the center of the distribution 
from which displacements are drawn by the known free- 
particle drift term. The translation into elastic collisions 
as in the flow-free case still ensures crgodicity of the algo- 
rithm and can be expected to be a reasonable approxima- 
tion for small jAt. Lees-Edwards boundary conditions 
allow to match the resulting linear velocity profile with 
the periodic images of the simulation box. 

For the simulations presented here, we chose a time 
step At = 0.01, resulting in Dq = 0.005. Initial config- 
urations have been allowed to equilibrate during runs of 
up to D t/a 2 = 4 x 10 4 , equivalent to 2 x 10 8 Brownian 
time steps. After equilibration, shear flow was instanta- 
neously switched on, and correlation functions have been 
measured for several waiting times t w thereafter. To im- 
prove statistics, the procedure has been repeated for 300 
independent runs at each density and shear rate. 

Figure [5] shows the self-intermediate scattering func- 
tions for various wave vectors in the gradient direction, 
at packing fraction tp = 0.79, and for a fixed shear rate 
7 = 0.02D /cr 2 (Pe = 0.02). The quiescent correla- 
tion functions are shown for comparison as dashed lines; 
they decay about two orders of magnitude slower than 
the shear-decorrelated ones, hence Pe ps 10 2 . Differ- 
ent waiting times are shown in dimensionless units jt w , 
the scaling expected from our theoretical observations 
above. One recognizes from the figure the same quali- 
tative trends as found in the Newtonian-dynamics sim- 
ulation of Ref. [l[ and consistent with our theory: the 
difference between the various waiting-time dependent 
correlation functions is most pronounced at intermedi- 
ate times, when the functions start to decay from their 
respective plateaus. Following switch-on, the transient 
correlation function stays close to its equilibrium counter- 
part up to jr ss O(0. 01), although deviations set in ear- 
lier for higher q. This is consistent with the picture that 
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FIG. 6: Tagged-particle density correlation functions 
from stochastic-dynamics computer simulation of a two- 
dimensional hard-disk system at area fraction tp = 0.79, just 
below the glass transition, and shear rate Peo = 0.02. Curves 
from top to bottom correspond to wave numbers qa = 1.5, 
4.8, 9.7, and 14.6, for wave vectors in the direction perpen- 
dicular to the flow direction. Thick solid and dashed lines are 
transient (t w = 0) and stationary (t w — ► oo) correlators, re- 
spectively. Thin lines represent different t w > as indicated, 
the dotted curves are quiescent equilibrium correlators. 

for fluctuations probing smaller length scales, smaller ac- 
cumulated strains arc needed to deviate from the qui- 
escent state. Increasing t w , the stationary correlator is 
approached for jt w > 0.1, somewhat later than in the 
schematic model discussed above. This results since in 
the simulation, also the startup stress u(t w ) approaches 
its steady-state value later than in the schematic model 
(see the discussion of Fig. [5]). The stationary correla- 
tion function deviates earlier from the equilibrium one 
than the transient one, and decays slower; generally all 
im-dependent correlation functions are found to merge 
again at longer times. We note in passing that the long- 
time decay of the transient correlation function could be 
fitted with a "compressed" exponential function, result- 
ing in exponents [i = 1.1, 1.8, 1.8, and 1.4 for the four 
different wave-vectors shown. In the corresponding MD 
simulation of Ref. [l[, values of fx ranging from 1.2 to 2.4 
have been found, increasing with increasing wave number 
P. 

Figure [7] displays the mean-squared displacements ob- 
tained from the ED-BD simulation, for fixed packing frac- 
tion tp = 0.79, but various shear rates 7 covering the 
regime Pe > 1 . Again in qualitative agreement with pre- 
vious Newtonian-dynamics results, and also with MSD 
curves obtained from confocal microscopy on colloidal 
suspensions [l|, for fixed 7, the curves for different t w 
all collapse for short and for long times onto the steady- 
state curve, deviating at intermediate times; the tran- 
sient MSD in the simulation deviates from the equilib- 
rium curve at jr ~ 0.02 and crosses over to the steady- 
state curve via a super-diffusive regime. Determining 
an effective exponent via the logarithmic derivative of 
the MSD, G?logt5r 2 (r)/<iiogT, for the largest 7 shown in 
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decreasing Ct w (r) for increasing t w shown in Fig. [2j i.e. 



FIG. 7: Mean-squared displacements in the shear-gradient 
direction for the stochastic-dynamics computer simulation of 
a hard disk system at area fraction (p — 0.79, as in Fig. [6] but 
for different Peclet numbers, Peo = 2 x 1CT 2 , 2 x 1CT 3 and 
2 x 10 -4 from left to right. Dotted line: quiescent equilibrium 
MSD. 
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FIG. 8: Waiting-time dependent tagged-particle density cor- 
relation functions for a two-dimensional hard-disk system un- 
dergoing stochastic dynamics, as in Fig.[6]for qa = 4.8, but for 
different waiting times corresponding to accumulated strains 
■yt w where a stress overshoot is seen, marked in the inset. 



Fig. [7] yields 8r 2 (r) r 1,9 , comparable with the expo- 
nent found in Ref. [Ij for the colloidal suspension, slightly 
smaller than the one extracted from the Newtonian dy- 
namics simulation (w 2.1). 

The results shown in Figs. [6] and [7] are to be compared 
to the schematic-model results shown in Figs. [2] and [5] 
Regarding the cross-over from transient to stationary cor- 
relation functions, the agreement is indeed qualitative. 
This holds despite the fact that, as mentioned above, the 
schematic model we employed misses the superdiffusivc 
regime in the MSD for technical reasons. 

In the schematic model, a(t w ) > is a monotonically 
increasing function of t w . We hence get a sequence of 



C twi (r)>C tw2 (r) 



for t Wl ^ t"W2,i 



(32) 



or reverse for the mean-squared displacement. This or- 
dering rule is well obeyed by the curves shown in Fig. [G] 
and the related ordering of the MSD is verified in Fig. [7] 
Note however that in simulations as well as in the mi- 
croscopic MCT the startup stress exhibits an overshoot 
connected to a small negative dip in the transient dynam- 
ical shear modulus []],[l3]], rendering a(t w ) nonmonotonic 
as a function of t w . Hence the ordering of the correlation 
functions given in Eq. (|32[) could in principle be violated 
for a small t^-window, compare Eq. (|23[) . In particu- 
lar, Ct w (r) for some fixed t w = 0(0.1/7) could conceiv- 
ably be smaller than the steady-state correlator. In the 
schematic model, this effect is not contained. 

However, as shown by Fig. [8l our simulation indeed 
indicates such a crossing of correlators as a function of 
t w . Here, waiting times were chosen to sample the 7^ 
region around the stress overshoot. For t w correspond- 
ing to the maximum stress, Ct m (r) is found to fall be- 
low the stationary correlator, while this is not the case 
for smaller t w . While lending credibility to our approxi- 
mations, Fig. [8] also indicates that the approximation in 
Eq. (|23|) of factorizing a(t w ) is an oversimplification. For 
example at jt w sa 0.04, the startup stress measured in 
the simulation reaches its steady-state value for the first 
time, before entering the overshoot region, and still, the 
corresponding correlator C tiu (r) differs from the station- 
ary Coo(r). 

The absence of the stress overshoot in our schematic 
model also implies that the stationary stress is reached 
earlier; waiting times of the order jt w sa 0.1 were suf- 
ficient to enter the stationary regime in Fig. [2] In the 
simulation, the corresponding t w are slightly larger, since 
one has to wait for the stress-overshoot region to be sur- 
passed. The latter causes the transient correlation func- 
tions to approach the steady-state ones only for jt w « 1 . 

Motivated by the qualitative agreement, we now turn 
to a more general test of Eq. (|23j) . by checking the pre- 
dicted relation among the three correlation functions 
(equilibrium, stationary, and finite t w ) for the simulation 
data. This will be done both for the Brownian-dynamics 
data set just discussed, and also for the Newtonian- 
dynamics simulation data found in Ref. [l|, in order to 
test the generality of our approximation regarding differ- 
ent forms of the short-time dynamics. 

In performing the comparison to follow, we are bur- 
dened by the fact that the calculation of Ct ro (r) from 
Cq(t) and C e (r) performed with experimental or simula- 
tion data is quite unstable, due to the roughening effect of 
the numerical derivative and a cancellation of small terms 
when all correlators are close to their plateau values. 
In the schematic model, Cq(t) and C e (r) were available 
with high enough precision. We therefore have to turn 
toward Eq. $%o§ and, for the MSD, Eq. (g7]). Although 
mathematically identical, these forms give fewer numer- 
ical difficulties, as the differentiation can be replaced by 
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a much smoother numerical integration. Unfortunately, 
this makes the procedure somewhat less intuitive: given 
C e (r) and Ct w (r) for some t w , we can now calculate the 
transient correlator Co(r) but will, due to the nature of 
the approximations involved or due to numerical inaccu- 
racies, get differing predictions for Co(t) from different 
t w . This will thus test the accuracy of our approxima- 
tions for different t w . 

In addition, a(t w ) appearing in Eqs. ([26]) and (|27p is 
in general not known or easily determined. However, for 
small t w , it can be replaced by its first order expansion, 
and we have a parameter-free prediction for Co (r) out of 
C e (r) and C tra (r) that we test in the next section. After 
that, we turn to the case t w — > oo, to demonstrate with 
cr(oo) taken as a (wave- vector independent) fit parameter 
the qualitative agreement of our result with the data. 



A. Small Waiting Times 

Figure [5] demonstrates the quality of our approxima- 
tion, Eq. (|2"3")l . for small waiting times, where o{t w ) = 
t w +0(t^,) holds. Both Brownian dynamics and MD sim- 
ulation results were used to probe the relation among the 
three correlator types via Eq. (f2"rJ)) . For the ED-BD simu- 
lation, the density is given by <p = 0.79, and Peo = 0.02, 
as in Fig. [51 For the MD, the temperature is given by 
T = 0.14 (close to the glass transition on the fluid side) 
and jy/ (m<j 2 )/e — 6 x 10 -4 , where e sets the particle 
interaction strength, see Ref. [1| for details. Note that in 
the figure, several wave vectors q in the gradient direc- 
tion were chosen. Recall that in Eqs. ([2lt]) and (j2"o]l , the 
wave-vector dependence is only implicit through that of 
C e (r) and C tvj {T). Since the projection onto stresses be- 
comes exact for small waiting times, see Eq. (fTT)) . Fig. [S] 
shows explicitely the accuracy of our approximation for 
the waiting time derivative, Eq. (|2"Tj) . Noting that this 
approximation was central for the derivation of nontrivial 
fluctuation-dissipation-ratios in Refs. [IJ (HI gives strong 
support for the validity of the results found there: Fig. [9] 
points out that this approximation captures well the ob- 
servable dependence, at least for the tagged-particle cor- 
relation functions we study. Let us emphasize again, that 
for this comparison, no free parameter appears in the 
equations, as a(t w ) = 1 holds for tagged-particle density 
fluctuations, and a(t w ) can be expanded to first order 
in t w . Nevertheless, both ED-BD and MD simulation 
data for the transient correlation function show excellent 
agreement with the curve calculated via Eq. ([26)) . 

A similarly good agreement is found for the mean- 
squared displacements in the directions perpendicular to 
the shear flow, as demonstrated by Fig. QUI again for both 
Brownian dynamics and MD simulation. In the latter, to 
improve statistics, the MSD has been averaged over both 
the y- and z-direction. The figure in particular high- 
lights that the failure to reproduce superdiffusive MSDs 
is completely within the schematic model we chose to il- 
lustrate the equations in Sec. IIV1 Similar to above, the 




log (t-t w ) y 

FIG. 9: (a) Parameter-free calculation of the transient cor- 
relation function from Eq. (|26[) from event-driven Brownian- 
dynamics (ED-BD) computer simulation of hard disks: dot- 
ted and dashed lines are the equilibrium and small-t^, corre- 
lation functions, respectively, at Pe « 100, Peo — 0.02 and 
jt-w = 0.010, with parameters as in Fig. [6] Solid lines in- 
dicate the resulting t w = transient as calculated from the 
approximation, Eq. ()26p . while symbols show the correspond- 
ing curves determined from computer simulation directly, (b) 
Same calculation, but for molecular-dynamics (MD) simula- 
tion data of Ref. pj, for a 3D binary soft-sphere mixture at 
Pe ~ 10 3 , -ytw — 0.017, and wave vector magnitudes qa = 2.3, 
6.0, and 12.3 (top to bottom) in the gradient direction. 



approximation for the waiting time derivative contained 
in Eq. (|25j) was central in deriving the nonequilibrium 
Einstein relation, Ref. [321 ] . The accuracy of this approx- 
imation as shown in Fig. 1101 hence supports the relations 
found there. 



B. Large waiting times 

Having verified that Eqs. (j23j) and ([25]) give a quan- 
titatively correct account for the relation between equi- 
librium, waiting-time dependent, and transient correla- 
tion function for short waiting times, we now turn to 
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FIG. 10: Parameter-free calculation of transient mean- 
squared displacements from equilibrium and sma,\\-t w data via 
Eq. (|27[) . parameters and symbols as in Fig. [9] (a) ED-BD 
computer simulation, (b) MD computer simulation. 



the t w — > oo limit, investigating the relationship be- 
tween stationary, equilibrium, and transient dynamics. 
We again follow the route outlined above, i.e., we de- 
termine the transient correlation functions and mean- 
squared displacements from the simulated equilibrium 
and stationary ones via Eqs. (|2l))) and (|2"T)) . and compare 
to the simulated transient correlation function in order 
to test the validity of our approach. 

To continue we then need cr(oo) entering our relation; 
this quantity could in principle be calculated following its 
definition in Eq. (fT4|) . In practice, we will use cr(oo) as a 
fit parameter, noting that it should be ^-independent (in 
general /-independent) to be meaningful. 

For the Brownian dynamics data (choosing the same 
state point and shear rate as in the small-i„, test above) , 
we find 7<r(oo) = 0.04 to give very satisfying results, 
shown in Fig. [IT] For the molecular-dynamics simula- 
tion data, 717(00) = 0.12 was fitted differently, as the 
parameter will depend not only on the shear rate but 
also the details of the interaction potential. Stationary 
correlators have been obtained from the simulation after 



FIG. 11: Calculation of the transient correlation function 
from equilibrium and steady state via Eq. (|26[) . Symbols and 
parameters are as in Fig. |9l but with the small-f^ curve re- 
placed by the t w — t 00 steady-state one. The parameter a(oo) 
from Eq. (|29[) has been fitted, (a) Comparison with ED-BD 
computer simulation data with 717(00) = 0.04. (b) Compari- 
son with MD computer simulation data with 7<r(oo) = 0.12. 



•yt w = 1, where no significant further change with t w was 
observed. As expected from the nature of our approxima- 
tion, which decouples two time dependences and is hence 
better for small t w , the agreement for the simulated and 
the calculated transient correlator is still very good, but 
somewhat less precise than in Fig. [9l Deviations in the 
comparison with Brownian dynamics simulation data are 
most pronounced for large q, where the predicted decay 
of the transient correlation function is too slow. 

For the molecular-dynamics data, larger deviations oc- 
cur, as exemplified in the lower panel of Fig. 111! The 
same qualitative trend with wave number holds as in the 
comparison with the ED-BD simulation, but the curve 
shape predicted for decay of the transient correlation 
function from the stationary one is markedly different 
from the one observed, especially in the initial decay 
from the plateau. Hence, while our approximation seems 
to hold for both ED-BD and MD equally well when ap- 
plied with small enough waiting times, as in Fig. [SJ it 
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FIG. 12: Calculation of transient mean-squared displace- 
ments in the gradient direction from equilibrium and station- 
ary curves via Eq. (|27[) with <r(oo)-values as in Fig. 1111 Curves 
correspond to the parameters as in Fig.[TTJ A dash-dotted line 
indicates a unit slope corresponding to diffusive motion, (a) 
ED-BD computer simulation, (b) MD computer simulation. 



is more adapted to ED-BD when applied with large t w . 
Let us point out that the time scale for the relaxation 
of the transient correlator is captured well in our ap- 
proximation although it differs appreciably from the one 
of the stationary correlator at the larger q shown. The 
difference between stationary and transient correlation 
functions in our theory is proportional to the parameter 
tr(oo). Changing it e.g. to smaller values gives tran- 
sient functions which are closer to the stationary ones in 
Figs. [TT1 (and also Fig. [12l see below). 

A similar picture arises from the MSD in the gradient 
direction perpendicular to the shear flow, shown again 
for both ED-BD and MD in Fig. [12] As in Fig. HU we 
have chosen jt w = 1 (ED-BD) and jt w tn 0.25 (MD) 
to approximate the stationary MSD in the simulations. 
The transient MSDs have been calculated using Eq. ([2"7]l , 
keeping the fit parameter ct(oo) fixed as determined from 
the comparison of the correlation functions above. The 
agreement is then seen to be qualitative again, with some 
quantitative deviations mostly regarding the logarithmic 



slope of the MSD at intermediate times. These deviations 
are larger for the MD data than for the ED-BD, and it 
will be interesting to explore further these differences in 
the waiting-time dependent evolution of sheared systems 
obeying different types of short-time dynamics. 

It appears that, while for the transient dynamics at 
small t w , possibly qualitative differences are found be- 
tween MD and BD (see the discussion above in conjunc- 
tion with Fig. [T0|) . this docs not seem to be the case for 
either the quiescent or the stationary case - the latter 
exemplified by the qualitative similarity of the large-t^, 
curves in Fig. 1121 



VI. DISCUSSION 

We have presented an approximate connection among 
the different dynamical two-point correlation functions 
characterizing the nonequilibrium dynamics following a 
sudden commencement of linear shear flow. It describes 
the evolution from the pre-shear equilibrium state to the 
stationary state attained under shear after sufficiently 
long waiting times. 

The approximations have been derived following the 
integration-through transients procedure as a general 
tool to evaluate averages involving the unknown nonequi- 
librium distribution function in terms of history integrals 
over the known equilibrium one. ITT is exact in prin- 
ciple, but in evaluating the two-time three-point corre- 
lation functions, we are forced to introduce approxima- 
tions. Here we have made the assumption that the dy- 
namical variables in whose evolution we are interested in, 
are, in their evolution from equilibrium to steady state, 
dominated by an overlap with the microscopic stress ten- 
sor, or more generally the fluctuations arising from ap- 
plying the nonequilibrium part of the time-evolution op- 
erator to the distribution function. 

An analysis of two different sets of computer simula- 
tion data close to the glass transition for various wave 
numbers, shear rates, and waiting times suggests that 
our approximative formulas do indeed capture the essen- 
tial features of the dynamical evolution from equilibrium 
to steady state. To test the generality of the discussed 
effects with respect to different short-time dynamics and 
dimensionality, we both analyzed previous 3D Newtonian 
molecular dynamics data and performed 2D event-driven 
Brownian dynamics (ED-BD) simulations leading to dif- 
fusive (overdamped stochastic) short-time motion. 

We have made plausible the approximations by refer- 
ring to a schematic model of mode-coupling theory that 
allows to calculate all the desired quantities numerically 
with high precision. The theoretical predictions are seen 
to broadly agree with previous results on both molecular- 
dynamics computer simulations and on colloidal experi- 
ments reported in Ref. [l[ , and also with the hard-sphere 
stochastic-dynamics simulations we performed. As a gen- 
eral trend, the transient correlation functions follow the 
equilibrium ones for times up to jr sa 0.01; after this 
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strain is reached, they deviate and decay much more 
rapidly in the regime of high dressed Peclet numbers we 
are interested in. For the tagged-particle density correla- 
tion functions, the deviation sets in earlier for larger wave 
number (tested in the gradient direction). For smaller Pe 
(or in the case of the MSD), the decisive strain is closer 
to 77 w 0.1: as mentioned before [l|, this value is sus- 
piciously close to a typical particle localization length, 
supporting the picture that the steep decay of the tran- 
sient correlation function marks the breaking of nearest- 
neighbor cages under shear. As the waiting time t w in- 
creases, this steep decay, found to be super-exponential 
in our model, broadens to settle on the t w — > oo station- 
ary limit for waiting times of yt w sa 1. For most cases, 
this gives rise to an ordering of ^-dependent correlators, 
with the transient, t w = 0, correlator being the largest 
(i.e. slowest), and the stationary, t w = oo, correlator be- 
ing the smallest (i.e. fastest) in an intermediate window 
of rescaled times jr. However, this ordering is violated 
for waiting times corresponding to strains jt w where a 
non-monotonic variation is observed in the shear stress 
after startup, the so called stress overshoot. Our approx- 
imation recovers this subtlety. 

As a side note, it should be mentioned that in the case 
of switching on steady shear, the waiting-time dependent 
correlation functions arc found to all decay on time scales 
of 0(1/7) when leaving the linear-response regime, i.e., 
when Pe > 1. This is in notable difference to waiting- 
time dependences usually discussed in the aging dynam- 
ics of glasses following a sudden temperature quench or 
in soft-matter systems [33U35| , where the final relaxation 
time strongly increases as a function of increasing t w , giv- 
ing rise to an ordering of correlators that is the reverse 
of what we discussed here. 

Let us remark that our schematic-model analysis sug- 
gests to discuss the peculiarities of the transient dynam- 
ics in colloidal and non-colloidal systems in separate as- 
pects: for both types of short-time motion "compressed 
exponential" relaxation is found in simulations, which in 
colloidal systems can be seen as the hallmark of nonequi- 
librium dynamics (as the equilibrium dynamics is con- 
fined to show at most exponential relaxation for any cor- 
relation function as long as only structural relaxation is 
involved). In fact, the noncquilibrium steady-state dy- 
namics according to our observations is again better de- 
scribed by pure exponentials, so that the compressed ex- 
ponentials may be the signature of noncquilibrium and 
non-stationary relaxation. This faster-than-exponential 
relaxation does, however, not directly translate into a 
super-diffusive (or even ballistic) regime in the mean- 
squared displacement, another feature observed in simu- 
lations on both types of system. As already investigated 
in detail in Ref . |l| , this latter feature can arise indepen- 
dently and is directly related to the stress overshoot. 

The mechanism by which the dynamics evolves from 
equilibrium via transient to stationary dynamics appears 
to have aspects that are independent on the details of 
the short-time motion: The points in time where de- 



viations are first seen from the equilibrium correlator, 
and last seen with respect to the steady state, are deter- 
mined by the total accumulated strain jr. The decisive 
value of 7T w 0.1 turns out to be suspiciously close to 
a typical cage size in supercooled liquids, as has been 
noted previously. We reiterate that this generality in 
the evolution from equilibrium to steady state does not 
mean that an equally strong statement holds as it does 
in equilibrium, where the long-time part of the dynamics 
itself is independent on the type of short-time dynam- 
ics. In fact, analyzing mean-squared displacements and 
tagged-particle density-correlation functions from both 
ED-BD and MD simulations, the MD ones appear to 
show a somewhat stronger transient effect than the ED- 
BD ones. Determining effective exponents for the MSD 
(from a logarithmic derivative) , the MD data gives larger 
deviations from sub-diffusive structural relaxation. The 
compressed-exponential exponents obtained from the in- 
coherent scattering function can be significantly larger 
than 2 in MD, but remain below 2 for the ED-BD data 
we analyzed. This occurs in a time window, where for 
the quiescent system structural relaxation prevails, and 
where one thus expects universal aspects of slow relax- 
ation that are independent of the short-time dynamics. 
It should be noted that the simulation models are dif- 
ferent, so that we cannot exclude the difference observed 
in the transients being due to details of the interparticle 
interactions. 

Our approximation correspondingly shows largest de- 
viations when relating equilibrium, transient, and sta- 
tionary correlation functions at large q for the molecular- 
dynamics simulation, even though it performs surpris- 
ingly well for the same MD data when restricting it to 
smaller waiting times. This may indicate that the pecu- 
liar form of the t w = transient in molecular dynamics 
is not fully captured by considering all dynamical relax- 
ations to be governed by the dynamics of the potential 
part of the local stress tensor (as we did to arrive at our 
approximation), but that non-potential parts of the dy- 
namics may play a role in this regime, too. 

Nevertheless, the overall quality of the waiting-time- 
derivativc approximation (Eq. (|21[) ) is found to be re- 
markably accurate. This lends further support to 
the earlier discussion of nonequilibrium fluctuation- 
dissipation ratios and Einstein relations within MCT and 
its schematic models, as the equivalent approximations 
have been used there [Tl], [ll| HH . 
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Appendix A: Hard Spheres 

For the hard sphere system, the instantaneous shear 
modulus diverges, rendering the projector V a leading to 
Eqs. (fl"3j) and (JT4J) ill-defined. For a proper treatment of 
this singular limit, one can use instead 

= e n1 S0 cx xy ){a xy e 2Q ' ' s ° <r xy )~ l {a xy e n] ' a ° (Al) 



for some fixed small sq. With this choice, 
Eq. (|13p remains unmodified, but the insertion of 
exp[f^ so] Va{so) exp[— fi' sq] changes Eq. (fT4|) to 



(A2) 



where the integrand remains finite for all t w > 2sq. The 
approximation thus, in the case of hard spheres, consists 
of coupling the variable / to the shear stress at some 
small cutoff time sq. 
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